import typing
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import ufl
import viskex
Generate and plot meshes
interval = dolfinx.mesh.create_unit_interval(mpi4py.MPI.COMM_WORLD, 6)
square_tria = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 3, 3, dolfinx.mesh.CellType.triangle)
square_quad = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 3, 3, dolfinx.mesh.CellType.quadrilateral)
cube_tetra = dolfinx.mesh.create_unit_cube(mpi4py.MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.tetrahedron)
cube_hexa = dolfinx.mesh.create_unit_cube(mpi4py.MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.hexahedron)
viskex.dolfinx.plot_mesh(interval)
viskex.dolfinx.plot_mesh(interval, dim=1)
viskex.dolfinx.plot_mesh(interval, dim=0)
viskex.dolfinx.plot_mesh(square_tria)
viskex.dolfinx.plot_mesh(square_tria, dim=2)
viskex.dolfinx.plot_mesh(square_tria, dim=1)
viskex.dolfinx.plot_mesh(square_tria, dim=0)
viskex.dolfinx.plot_mesh(square_quad)
viskex.dolfinx.plot_mesh(square_quad, dim=2)
viskex.dolfinx.plot_mesh(square_quad, dim=1)
viskex.dolfinx.plot_mesh(square_quad, dim=0)
viskex.dolfinx.plot_mesh(cube_tetra)
viskex.dolfinx.plot_mesh(cube_tetra, dim=3)
viskex.dolfinx.plot_mesh(cube_tetra, dim=2)
viskex.dolfinx.plot_mesh(cube_tetra, dim=1)
viskex.dolfinx.plot_mesh(cube_tetra, dim=0)
viskex.dolfinx.plot_mesh(cube_hexa)
viskex.dolfinx.plot_mesh(cube_hexa, dim=3)
viskex.dolfinx.plot_mesh(cube_hexa, dim=2)
viskex.dolfinx.plot_mesh(cube_hexa, dim=1)
viskex.dolfinx.plot_mesh(cube_hexa, dim=0)
Generate and plot subdomains
def mark_subdomains(mesh: dolfinx.mesh.Mesh) -> dolfinx.mesh.MeshTags:
"""Mark left and right subdomains in a given mesh with values 1 and 2, respectively."""
def left_subdomain(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Condition that defines the left subdomain."""
return x[0] <= 1.0 / 3.0 # type: ignore[no-any-return]
def right_subdomain(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Condition that defines the right subdomain."""
return x[0] >= 2.0 / 3.0 # type: ignore[no-any-return]
subdomains_entities = dict()
subdomains_values = dict()
for (subdomain, subdomain_id) in zip((left_subdomain, right_subdomain), (1, 2)):
subdomains_entities[subdomain_id] = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, subdomain)
subdomains_values[subdomain_id] = np.full(subdomains_entities[subdomain_id].shape, subdomain_id, dtype=np.int32)
subdomains_entities_unsorted = np.hstack(list(subdomains_entities.values()))
subdomains_values_unsorted = np.hstack(list(subdomains_values.values()))
subdomains_entities_argsort = np.argsort(subdomains_entities_unsorted)
subdomains_entities_sorted = subdomains_entities_unsorted[subdomains_entities_argsort]
subdomains_values_sorted = subdomains_values_unsorted[subdomains_entities_argsort]
subdomains = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, subdomains_entities_sorted, subdomains_values_sorted)
return subdomains
interval_subdomains = mark_subdomains(interval)
square_tria_subdomains = mark_subdomains(square_tria)
square_quad_subdomains = mark_subdomains(square_quad)
cube_tetra_subdomains = mark_subdomains(cube_tetra)
cube_hexa_subdomains = mark_subdomains(cube_hexa)
viskex.dolfinx.plot_mesh_entities(
interval, 1, "subdomains", interval_subdomains.indices[interval_subdomains.values == 2])
viskex.dolfinx.plot_mesh_entities(
interval, 1, "subdomains", interval_subdomains.indices, interval_subdomains.values)
viskex.dolfinx.plot_mesh_tags(interval, interval_subdomains, "subdomains")
viskex.dolfinx.plot_mesh_entities(
square_tria, 2, "subdomains", square_tria_subdomains.indices[square_tria_subdomains.values == 2])
viskex.dolfinx.plot_mesh_entities(
square_tria, 2, "subdomains", square_tria_subdomains.indices, square_tria_subdomains.values)
viskex.dolfinx.plot_mesh_tags(square_tria, square_tria_subdomains, "subdomains")
viskex.dolfinx.plot_mesh_tags(square_quad, square_quad_subdomains, "subdomains")
viskex.dolfinx.plot_mesh_tags(cube_tetra, cube_tetra_subdomains, "subdomains")
viskex.dolfinx.plot_mesh_tags(cube_hexa, cube_hexa_subdomains, "subdomains")
Generate and plot boundaries
def mark_boundaries(mesh: dolfinx.mesh.Mesh, subdomains: dolfinx.mesh.MeshTags) -> dolfinx.mesh.MeshTags:
"""
Mark internal and boundary facets in a given mesh with four different values.
Internal facets of left and right subdomains are associated with values 1 and 2, respectively.
Furthermore, boundary facets on the left and right boundaries are associated with values 3 and 4,
respectively.
"""
def left(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Condition that defines the left boundary."""
return abs(x[0] - 0.) < np.finfo(float).eps # type: ignore[no-any-return]
def right(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Condition that defines the right boundary."""
return abs(x[0] - 1.) < np.finfo(float).eps # type: ignore[no-any-return]
boundaries_entities = dict()
boundaries_values = dict()
for (boundary, boundary_id) in zip((left, right), (3, 4)):
boundaries_entities[boundary_id] = dolfinx.mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, boundary)
boundaries_values[boundary_id] = np.full(
boundaries_entities[boundary_id].shape, boundary_id, dtype=np.int32)
boundaries_entities_unsorted = np.hstack(list(boundaries_entities.values()))
boundaries_values_unsorted = np.hstack(list(boundaries_values.values()))
boundaries_entities_argsort = np.argsort(boundaries_entities_unsorted)
boundaries_entities_sorted = boundaries_entities_unsorted[boundaries_entities_argsort]
boundaries_values_sorted = boundaries_values_unsorted[boundaries_entities_argsort]
boundaries = dolfinx.mesh.meshtags(
mesh, mesh.topology.dim - 1,
boundaries_entities_sorted, boundaries_values_sorted)
cell_to_facets_connectivity = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim - 1)
new_facets_indices = list()
new_facets_values = list()
for cell in range(cell_to_facets_connectivity.num_nodes):
facets = cell_to_facets_connectivity.links(cell)
subdomain_index = np.where(subdomains.indices == cell)[0]
if subdomain_index.size > 0:
assert subdomain_index.size == 1
subdomain_value = subdomains.values[subdomain_index[0]]
if subdomain_value in (1, 2):
for facet in facets:
if facet not in boundaries.indices and facet not in new_facets_indices:
new_facets_indices.append(facet)
new_facets_values.append(subdomain_value)
boundaries_and_interfaces_entities_unsorted = np.hstack((boundaries.indices, new_facets_indices)).astype(np.int32)
boundaries_and_interfaces_values_unsorted = np.hstack((boundaries.values, new_facets_values)).astype(np.int32)
boundaries_and_interfaces_entities_argsort = np.argsort(boundaries_and_interfaces_entities_unsorted)
boundaries_and_interfaces_entities_sorted = boundaries_and_interfaces_entities_unsorted[
boundaries_and_interfaces_entities_argsort]
boundaries_and_interfaces_values_sorted = boundaries_and_interfaces_values_unsorted[
boundaries_and_interfaces_entities_argsort]
boundaries_and_interfaces = dolfinx.mesh.meshtags(
mesh, mesh.topology.dim - 1,
boundaries_and_interfaces_entities_sorted, boundaries_and_interfaces_values_sorted)
return boundaries_and_interfaces
interval_boundaries = mark_boundaries(interval, interval_subdomains)
square_tria_boundaries = mark_boundaries(square_tria, square_tria_subdomains)
square_quad_boundaries = mark_boundaries(square_quad, square_quad_subdomains)
cube_tetra_boundaries = mark_boundaries(cube_tetra, cube_tetra_subdomains)
cube_hexa_boundaries = mark_boundaries(cube_hexa, cube_hexa_subdomains)
viskex.dolfinx.plot_mesh_entities(
interval, 0, "boundaries", interval_boundaries.indices[interval_boundaries.values == 2])
viskex.dolfinx.plot_mesh_entities(
interval, 0, "boundaries", interval_boundaries.indices, interval_boundaries.values)
viskex.dolfinx.plot_mesh_tags(interval, interval_boundaries, "boundaries")
viskex.dolfinx.plot_mesh_entities(
square_tria, 1, "boundaries", square_tria_boundaries.indices[square_tria_boundaries.values == 2])
viskex.dolfinx.plot_mesh_entities(
square_tria, 1, "boundaries", square_tria_boundaries.indices, square_tria_boundaries.values)
viskex.dolfinx.plot_mesh_tags(square_tria, square_tria_boundaries, "boundaries")
viskex.dolfinx.plot_mesh_tags(square_quad, square_quad_boundaries, "boundaries")
viskex.dolfinx.plot_mesh_tags(cube_tetra, cube_tetra_boundaries, "boundaries")
viskex.dolfinx.plot_mesh_tags(cube_hexa, cube_hexa_boundaries, "boundaries")
Interpolate and plot scalar functions
def prepare_scalar_field_cases( # type: ignore[no-any-unimported]
mesh: dolfinx.mesh.Mesh,
expression: typing.Callable[
[typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]],
typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]
]
) -> typing.Tuple[dolfinx.fem.Function, typing.Tuple[ufl.core.expr.Expr, dolfinx.fem.FunctionSpace]]:
"""Prepare scalar field cases."""
scalar_function_space = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))
scalar_field = dolfinx.fem.Function(scalar_function_space)
scalar_field.interpolate(expression)
scalar_field_ufl = expression(ufl.SpatialCoordinate(mesh))
return scalar_field, (scalar_field_ufl, scalar_function_space)
interval_scalar_field, interval_scalar_field_ufl = prepare_scalar_field_cases(
interval, lambda x: x[0]**3)
square_tria_scalar_field, square_tria_scalar_field_ufl = prepare_scalar_field_cases(
square_tria, lambda x: x[0]**3 + x[1]**2)
square_quad_scalar_field, square_quad_scalar_field_ufl = prepare_scalar_field_cases(
square_quad, lambda x: x[0]**3 + x[1]**2)
cube_tetra_scalar_field, cube_tetra_scalar_field_ufl = prepare_scalar_field_cases(
cube_tetra, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
cube_hexa_scalar_field, cube_hexa_scalar_field_ufl = prepare_scalar_field_cases(
cube_hexa, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
viskex.dolfinx.plot_scalar_field(interval_scalar_field, "scalar")
viskex.dolfinx.plot_scalar_field(interval_scalar_field_ufl, "scalar")
viskex.dolfinx.plot_scalar_field(square_tria_scalar_field, "scalar")
viskex.dolfinx.plot_scalar_field(square_tria_scalar_field, "scalar", warp_factor=1.0)
viskex.dolfinx.plot_scalar_field(square_tria_scalar_field_ufl, "scalar")
viskex.dolfinx.plot_scalar_field(square_quad_scalar_field, "scalar")
viskex.dolfinx.plot_scalar_field(square_quad_scalar_field, "scalar", warp_factor=1.0)
viskex.dolfinx.plot_scalar_field(square_quad_scalar_field_ufl, "scalar")
viskex.dolfinx.plot_scalar_field(cube_tetra_scalar_field, "scalar")
viskex.dolfinx.plot_scalar_field(cube_tetra_scalar_field_ufl, "scalar")
viskex.dolfinx.plot_scalar_field(cube_hexa_scalar_field, "scalar")
viskex.dolfinx.plot_scalar_field(cube_hexa_scalar_field_ufl, "scalar")
Interpolate and plot vector functions
def prepare_vector_field_cases( # type: ignore[no-any-unimported]
mesh: dolfinx.mesh.Mesh,
expression: typing.Callable[
[typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]],
typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]
]
) -> typing.Tuple[dolfinx.fem.Function, typing.Tuple[ufl.core.expr.Expr, dolfinx.fem.FunctionSpace]]:
"""Prepare vector field cases."""
vector_function_space = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2))
vector_field = dolfinx.fem.Function(vector_function_space)
vector_field.interpolate(expression)
vector_field_ufl = ufl.as_vector(expression(ufl.SpatialCoordinate(mesh)))
return vector_field, (vector_field_ufl, vector_function_space)
square_tria_vector_field, square_tria_vector_field_ufl = prepare_vector_field_cases(
square_tria, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
square_quad_vector_field, square_quad_vector_field_ufl = prepare_vector_field_cases(
square_quad, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
cube_tetra_vector_field, cube_tetra_vector_field_ufl = prepare_vector_field_cases(
cube_tetra, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
cube_hexa_vector_field, cube_hexa_vector_field_ufl = prepare_vector_field_cases(
cube_hexa, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector")
viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector", glyph_factor=0.1)
viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector", warp_factor=1.0)
viskex.dolfinx.plot_vector_field(square_tria_vector_field_ufl, "vector")
viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector")
viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)
viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector", warp_factor=1.0)
viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector")
viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector", glyph_factor=0.1)
viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector", warp_factor=1.0)
viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector")
viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector", glyph_factor=0.1)
viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector", warp_factor=1.0)